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In hierarchical cosmological models^, galaxies grow in mass through the continual accretion 
of smaller ones. The tidal disruption of these systems is expected to result in loosely bound 
stars surrounding the galaxy, at distances that reach 10 — 100 times the radius of the cen- 
tral disk^'^. The number, luminosity and morphology of the relics of this process provide 
significant clues to galaxy formation history"', but obtaining a comprehensive survey of these 
components is difficult because of their intrinsic faintness and vast extent. Here we report a 
panoramic survey of the Andromeda galaxy (M31). We detect stars and coherent structures 
that are almost certainly remnants of dwarf galaxies destroyed by the tidal field of M31. An 
improved census of their surviving counterparts implies that three-quarters of M31's satel- 
lites brighter than My < —6 await discovery. The brightest companion, Triangulum (M33), 
is surrounded by a stellar structure that provides persuasive evidence for a recent encounter 
with M31. This panorama of galaxy structure directly confirms the basic tenets of the hi- 
erarchical galaxy formation model and reveals the shared history of M31 and M33 in the 
unceasing build-up of galaxies. 

Precise measurements of stars and star clusters within the Milky Way have contributed sig- 
nificantly to the development of a cosmological model for galaxy formation^'^. The discovery 
that the Sagittarius dwarf galaxy was being cannibalized by the Milky Way'' brought into sharp 
focus the role of satellite accretion in the build-up of a galaxy's mass. Models now propose that 
galaxies form within dark-matter haloes that grow through the continual accretion and merger of 
smaller sub-haloes. The predicted number of sub-haloes is at least a few orders of magnitude more 
than the number of dwarf galaxies observed as satellites around the Milky Way®'^. If we assume 
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that the underlying cosmology is correct, this implies that either significant numbers of satellites 
remain undiscovered or only a fraction of sub-haloes contain baryons (stars and gas). In either 
case, the number, luminosity and spatial distributions of satellite galaxies are important but poorly 
determined quantities whose values depend strongly on the processes through which baryons are 
accreted and retained by sub-haloes^°'^^. Many of the luminous sub-haloes are expected to be per- 
turbed, and even shredded, by the tidal field of the host galaxy, leaving behind stellar debris in the 
form of streams and substructures within a diffuse stellar halo^'^. 

Systematic studies of the Milky Way's stellar halo, such as with the Sloan Digital Sky 
Survey^^, have recently revealed a large number of dwarf galaxies and tidal streams. However, 
our viewpoint from within the Milky Way introduces selection effects and difficulties in interpre- 
tation, making a homogeneous global study difficult. To this end, we have initiated the "Pan- 
Andromeda Archaeological Survey" (PAndAS), a programme using the 1 -square-degree field-of- 
view MegaPrime/MegaCam camera on the 3.6 m Canada-France-Hawaii Telescope (CFHT). We 
are imaging the closest spiral galaxy, M31, and its less massive companion M33. Once completed, 
the survey will cover more than 300 square degrees (more than 70,000 kpc^) and extend to a max- 
imum projected radius from M3rs centre of Vp < 150 kpc. It is the largest contiguous imaging 
survey of a massive galaxy and spans the stellar halo out to extremely large radii. PAndAS builds 
on earlier Isaac Newton Telescope and CFHT photometric studies of this galaxy^^"^^. It surveys 
in the g (4, 140 - 5, 600A) and i (7, 020 - 8, 530 A) bands and resolves individual stars in M31 to 
depths of 51 = 26.5, i = 25.5 at a signal-to-noise ratio of 10. The programme started in August 
2008 and will continue until January 201 1. 
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Figure 1 shows the spatial density distribution of sources in our extant PAndAS fields that 
are consistent with red giant branch (RGB) stars in M31. Although the region surveyed so far 
covers ~ 220 square degrees (nearly 100-fold the area of the classical optical disk of M31), we 
find RGB stars everywhere across our survey. These stars trace the low-luminosity structure of the 
galaxy and reveal the vast extent of M31, challenging the commonly held impression of the size 
of this typical bright galaxy (and, by extension, other galaxies of similar luminosities). These stars 
are unlikely to have formed in situ at these radii because it is improbable that the local density of 
gas was high enough to promote star formation. Instead, it is expected that the stars have been 
accreted from dwarf galaxies or proto-galactic fragments, a conclusion that is consistent with the 
basic tenets of hierarchical galaxy formation. 

Our interpretation of halo stars as accreted relics is supported by the presence of multiple, 
large, coherent substructures over our survey area. Some of these features, such as the giant stellar 
stream^^ (no. 5 in Fig. 1), were previously known. New structures discovered in our survey 
include a radial overdensity along the northwest minor axis extending nearly 100 kpc from M31 
(no. 6), a diffuse structure to the southwest coherent over an arc spanning ~ 40 kpc at a distance 
of ~ 100 kpc from M31 (no. 7), and an apparent continuation of a previously known stream in the 
east^^, ~ 50 kpc from M31 that loops around to the north of the galaxy (no. 4). The large scale 
of the new structures is striking, as is their distance from the centre of M31. They are expected 
to maintain coherence for at least a few gigayears in the inner halo (r^ < 50 kpc), and at least a 
Hubble time in the outer halo (r^ > 100 kpc)^^. The implication is that these structures are the 
remains of previously accreted dwarf galaxies. None of the newly discovered stellar structures 
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clearly correlate with HI detections in the environs of M31 (ref. 19). Although the lack of young 
stars tracing sites of star formation implies that high HI column densities are not expected, more 
detailed predictions on the HI content of these substructures are not yet possible. 

The surviving counterparts of the progenitors of these substructures - the dwarf galaxies - 
are visible in Fig. 1 as concentrated, round, overdensities, and the large area surveyed permits 
an improved census of these objects. In Fig. 2a we plot the projected radial number density 
distribution of galaxies around M31; the profile shows no sign of declining within 150 kpc. In Fig. 
2b we plot the cumulative luminosity distribution of all galaxies around M31. We find that it is 
well described by the usual Schechter function^'^ with a faint-end slope of a = 0.98 ± 0.07, where 
the uncertainty represents the standard deviation. These distributions imply that 6 ± 4 satellites 
more luminous than My ~ —6 still remain to be discovered within 150 kpc (lOdegrees) of M31. 
Extrapolating the observed flat number-density profile outwards then suggests that M3 1 may have 
as many as 88 ± 20 such satellites out to Vp < 300 kpc (roughly equal to the expected virial radius 
of M3rs dark-matter halo; see Supplementary Information), only about one-quarter of which are 
currently known. The derived number is an upper limit, because we expect the radial profile to 
decline beyond the survey area. Thus, even accounting for observational incompleteness, there is 
still an order of magnitude too few satellites compared with the expected number of dark-matter 
haloes. If some of these haloes do not contain stars, then a comprehensive inventory of those that 
do may shed much light on the solution to the "missing satellites problem". 

The brightest of M3rs satellite companions, M33, is surrounded by a previously unknown 
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prominent stellar structure (no. 1 in Fig. 1). This feature has an extension stretching ~ 2 degrees 
(~ 30 kpc projected) to the northwest towards M3 1, nearly three times farther out than the classical 
disk of M33. A second extension is also visible in the south. One possible origin of this structure 
is that it is the accreted remains of a dwarf satellite, similar to the structures observed around M31. 
However, a long-standing puzzle about M33 is the existence of an extremely warped HI disk, with 
no apparent counterpart in the "pristine" stellar disk^^'^^. What process could distort the HI disk 
but leave the stellar disk unaffected? The newly discovered feature has a similar orientation to 
the HI warp, suggesting that we have discovered its optical counterpart. Furthermore, the general 
northwest - southeast symmetry of the distortion, and the broad alignment of the northwest exten- 
sion with the measured transverse velocity of M33 (ref. 23), are evidence of a tidal disturbance 
excited as this galaxy orbits around M31. 

We test this hypothesis by considering the M33 - M31 orbit. If the orbit of M33 carries it 
close to M3 1 , the tidal field of the latter will excite the M33 disk and eject stars, but the orbit cannot 
be so close that M33 is severely distorted or disrupted^^. We use the techniques described in the 
Supplementary Information to search for M33 orbits that are consistent with the known constraints 
and have reasonably close encounters with M31. The search reveals that the smallest possible 
pericentre distance is ~ 40 kpc, with close passage occurring a few billion years ago. We select 
several representative orbits and carry out n-body simulations^^ of the encounter. 

Figure 3 and Supplementary Movie 1 show results from one plausible interaction model. 
This orbit reproduces with good accuracy the observed distances^^, angular positions and radial 
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velocities of M31 and M33, as well as the proper motion of M33 (ref. 23). The interaction excites 
tidal tails in the M33 disk that wind up to form an extended distribution of stars with dimensions 
similar to those of the debris observed in the PAndAS data. In addition, the debris is warped away 
from the disk in the same sense as the known gaseous warp^^'^^, lending further support to the 

interaction hypothesis. 

Finally, we note that the encounter between M31 and M33 leaves an imprint on the larger 
system, exciting a mild warp and disturbances in the M3 1 disk at large radii. These phenomena are 
consistent with some of the unusual features observed in the outer regions of M31 (ref. 14), par- 
ticularly with the discovery of younger stellar populations in these fields^^ and with the measured 
rotational signature shared by many of the substructures^*^. It is plausible that these stars were 
originally formed in the thin disk and excited to their present locations by a galactic interaction, 
perhaps with M33 at an earlier phase in its orbit. The unrivalled panorama of galaxy structure pre- 
sented here for M31 therefore reveals the continuing role of accretion and interactions in shaping 
its properties, and is a startling visual demonstration of the truly vast scale of galaxies. 



8 




9 




10 




11 



Supplementary Information 

1. Galaxy Models 

We construct realistic, self -consistent models for M31 and M33 for use in our simulations. 
They consist of an exponential disk, a Sersic bulge, and a spherical cuspy dark matter halo. From 
these we construct a self-consistent distribution function and so generate equilibrium N-body sys- 
tems for simulations. The key structural parameters for the model used in the simulation are given 
in Supplementary Table 1. The masses of the three components, Mh, M^, and Mb are given in 
units of 10^° Mq. M33 is naturally truncated by the tidal field of M31 at a radius of ~50 kpc at 
the initial distance of 230 kpc in the model so we quote the M33 halo mass within this tidal radius. 
Note that for an isolated galaxy with M33's rotation curve, we would normally infer a virial radius 
> 100 kpc and virial mass about twice that quoted in the table below so it is important to truncate 
the model in this context. The M3 1 model extends to a virial radius of Vmriai — 280 kpc and we 
quote the halo mass inside this radius. The halo scale radius, r^, the disk scale length, R^, and the 
effective radius of the bulge. Re, are given in kiloparsecs. The Sersic index of the bulge is Ug. We 
have chosen a suitable mass-to-light ratio for the disk components to ensure that the disks are both 
intrinsically stable against the formation of a strong bar. To probe the evolution of the outer disks, 
we truncate the M33 model at 10 scale-lengths (13 kpc) and the M31 model at 7 scale-lengths (40 
kpc). These radii are larger than the classical size of the optical disks that extend to approximately 
5 scale-lengths. 

2. Orbital Model 
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We set our galaxy models up on a orbit that will lead to an interaction while satisfying 
the observed constraints. We define the orbit relative to the M31 centre. The assumed distances 
imply a current separation of M31 and M33 of 206 kpc. We move from heliocentric coordinates 
to a frame with an origin at the Galactic centre and correct for the effect of Galactic rotation 
using standard values of the rotation velocity and motion with respect to the Local Standard of 
Rest. The proper motion of M31 is as yet unmeasured so the transverse velocity components of 
M3 1 are unconstrained. However, if we assume values of the M3 1 proper motion then we can 
fully determine the current relative position and velocity and so determine the orbit. Finally, the 
inclination and position angles and direction of rotation of M3 1 and M33 determine the orientation 
of their respective disks in space. We assume that relative orientations of the two galaxies remain 
essentially unchanged during an interaction implying that the initial orientation of the two disks is 
the same as it is at the current epoch in an orbital model. 

We now try to find an orbit (or family of orbits) that satisfy the observed constraints. The 
main unknowns are the M31 transverse velocity components so we allow these to vary as free 
parameters. We first search the orbital space by integrating the motion of a M33 test particle in 
the gravitational potential of the M31 model defined above. Our goal is to find orbits with initial 
conditions that permit a strong enough tidal encounter to excite the observed features in M33. We 
therefore search for orbits that have close encounters. In test particle integrations, the observed 
constraints appear to limit the pericenter distance to a value no smaller than Vp ^ 40 kpc. In an 
effort to excite the strongest response we therefore confined ourselves to orbits with rp < 50 kpc 
which narrows the range of possible orbits considerably. Most orbits tend to be confined to a plane 
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nearly parallel to our line of sight and most of the interactions occur with the past few billion years. 
The model presented in this paper is representative of this family. 

Dynamical friction is a further complication. Its effect was significant and affected our pre- 
dictions of initial conditions based on the test particle integrations above. We therefore incorporate 
dynamical friction into our orbital search calibrating the parameters in the dynamical friction for- 
mula using the N-body simulations. In this way, it was possible to predict initial conditions for the 
simulations that would result in the M31 and M33 systems ending up in the current positions and 
relative velocities within the errors of the known constraints and assumed proper motion of M31. 

3. N-body Model 

We construct N-body models with a total of 6.7 million particles divided up according to 
Supplementary Table 2. This numerical resolution was adequate to illustrate the main features of 
a tidal encounter though future work will increase particle numbers to resolve finer details of the 
interaction. The simulation was carried out using a parallelized tree-code. 

4. Uniqueness of the model 

The observed constraints on the relative position and velocity of M33 along with the as- 
sumption of a close encounter lead to a range of models similar to the one presented in this paper. 
Another implicit assumption in determination of this orbit is the gravitational potential of M3 1 
with parameters described above. The outer mass profile and virial mass of M31 might be differ- 
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ent and this would change the detailed character of the orbit and interaction. However, the inner 
mass profile of M3 1 is well-constrained so the strength of the tidal encounter will not be signifi- 
cantly different in models where the outer mass profile varies. So while the model presented here 
is not unique it demonstrates with reasonable accuracy the effect of a recent close encounter on the 
evolution of M33 in a plausible orbit. 

The model demonstrates that such an encounter will excite tidal tails in the outer disk of M33 
shortly after the encounter as well as a warp. As the galaxy proceeds along its orbit, these features 
wind up leaving a disturbed outer disk. The extent of these features will depend in detail on the 
initial truncation radius of the disk and pericentric separation. Also, the phase of disk extension 
and warp will depend on the exact timing of the orbit. Matching up all these details simultaneously 
with all constraints is a difficult challenge, but the model presented here demonstrates quite clearly 
that a disturbance in the stellar disk and a warp with similar properties to the observed features can 
be generated by a recent interaction that is sufficiently close. 
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Figure 1 Stellar density map of Andromeda - Triangulum. A tangent-plane projection 
of the density distribution of stellar sources in our extant PAndAS imaging is shown, with 
colours and magnitudes consistent with RGB stars at the distance of M31 . The inset 
shows the central parts of our survey at higher resolution. Dashed circles represent the 
maximum projected radii of 150 and 50 kpc from M31 and M33, respectively. Scale im- 
ages of the disks of M31 and M33 are overlaid. Visible dwarf satellites are indicated with 
roman numerals. Numbers in circles indicate the largest and most obvious substructures 
detected in a visual inspection of the image: 1 , M33 structure; 2, 125-kpc stream (stream 
A); 3, stream C; 4, eastern arc (stream D); 5, giant stellar stream; 6, northwest minor-axis 
stream; 7, southwest cloud. Features 1, 6 and 7 (and part of 4) are new discoveries. 
Stellar sources were identified by using star-galaxy classification techniques described 
previously^"^'^^. Candidate RGB stars were selected by their position in a colourmagnitude 
diagram relative to theoretical isochrones^^ for a 12-Gyr population at the distance of M31 
(ref. 26), using only stellar sources with io < 23.5. A projection of the stellar density dis- 
tribution within a putative metallicity range -2.5 < [Fe/H] < -0.6 dex was created with 
0.02 X 0.02 pixels, smoothed with a Gaussian filter with dispersion 3 arcmin, and displayed 
with square-root scaling. The inset was created by combining, with a red-green-blue 
colour scheme, three such maps of M31 with metallicities of -0.4 < [Fe/H] < +0.2 dex, 
-1.3 < [Fe/H] < -0.4 dex and -2.3 < [Fe/H] < -1.3 dex, respectively (each with 
0.01 X 0.01 degree pixels, smoothed with a Gaussian filter with 2 arcmin dispersion) and 
displayed with logarithmic scaling. Not all structures are visible with this (or any other) 
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choice of metallicity cut, filter and scaling. 

Figure 2 Distribution of M31 dwarf galaxies, a, Points show the projected radial number 
density of dwarfs (E), derived for galaxies within 10 degrees of M31 in the survey region. 
Horizontal error bars show the size of the bin, with points plotted at the mean galaxy 
radius in that bin. Vertical error bars are Poissonian. Galaxy names in black indicate the 
positions of the galaxies contributing to each bin; the remaining galaxies are named in 
grey. The radial ranges of the substructures indicated in Fig. 1 are also shown, b. The 
luminosity distribution, using all galaxies within the survey. Error bars are Poissonian. 
The dashed line shows the best-fit Schechter function, with slope a = 0.98 ± 0.07. Dark 
and light grey areas show ^a and 2a deviations in a, respectively. In both panels we are 
probably missing some satellites within ~ 4 degrees of M31 , where the high stellar density 
makes it difficult to detect faint satellites. We can detect galaxies brighter than My < 26, 
although derivation of the exact incompleteness levels remains to be calculated for the 
completed survey. To avoid these selection effects, the slope of the luminosity function is 
derived by using a fit of a Schechter function to galaxies brighter than My < -8, where 
incompleteness is not an issue, and we note that the extrapolation to lower luminosities 
fits these data well. Because of the lack of bright galaxies, we cannot independently 
derive M*. Instead, we fix M* at a range of values and derive corresponding values of a, 
finding that a is robust to within la for -19.5 > M* > -21.5 (with a trend such that fainter 
values for M* correspond to lower values for a). The quoted value of a corresponds to 
M* ^ -21. 
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Figure 3 An M31 - M33 interaction, a, The evolution of M33 about M31 along an orbit 
consistent with the angular positions, distances^^ and radial velocities of M31 and M33 
and with M33's proper motion^^. Here, M33 starts 3.4 Gyr ago at a distance of r < 200 kpc 
on the far side of M31 falling down the line of sight to the Milky Way. After 800 Myr, M33 
reaches pericentre and proceeds across our line of sight towards the southeast, reaching 
apocentre about 900 Myr ago before falling back towards M31 to its current position. 
Dots tracing the orbit are separated by 49 Myr to give a sense of the speed along the 
orbit. The lower inset shows a perpendicular view of M33's orbit, b. Quantification of 
the expected K-band surface brightness of M33 at different times from face-on and edge- 
on perspectives. The inner red/orange region with i^k < 22 magarcsec"^ defines the 
size of the usual optical disk of M33 seen in images. The initial equilibrium models^° for 
the two galaxies consist of a disk, bulge and dark-matter halo with structural parameters 
that accurately reproduce the observed surface brightness profiles and rotation curves. 
Because the mass profile of M31 beyond 100 kpc is not well constrained by observations, 
we appeal to cosmological arguments^"^ that predict a mass of 2.5 x 10^^ solar masses 
within r < 280 kpc. See the Supplementary Information for more details. 
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Halo 


Disk 

Md Rd 


Bulge 

Mb Re Us 


M33 
M31 


8 12 
247 22 


0.26 1.3 
5.8 5.8 


0.055 0.84 2.5 
3.3 1.5 1.0 



Table 1 : Equilibrium galaxy parameters 





Halo 


Disk 


Bulge 




Nh 


Nd 


Nb 


M33 


2M 


1M 


200K 


M31 


2M 


1M 


500K 



Table 2: Numerical resolution 
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